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ABSTRACT 

Numerical simulations of weakly magnetized and strongly magnetized rela- 
tivistic jets embedded in a weakly magnetized and strongly magnetized stationary 
or weakly relativistic {v = c/2) sheath have been performed. A magnetic field 
parallel to the flow is used in these simulations performed by the new GRMHD 
numerical code RAISHIN used in its RMHD configuration. In the numerical 
simulations the Lorentz factor 7 = 2.5 jet is precessed to break the initial equi- 
librium configuration. In the simulations sound speeds are < c/ \/3 in the weakly 
magnetized simulations and < 0.3 c in the strongly magnetized simulations. The 
Alfven wave speed is < 0.07 c in the weakly magnetized simulations and < 0.56 c 
in the strongly magnetized simulations. The results of the numerical simula- 
tions are compared to theoretical predictions from a normal mode analysis of 
the linearized relativistic magnetohydrodynamic (RMHD) equations capable of 
describing a uniform axially magnetized cylindrical relativistic jet embedded in 
a uniform axially magnetized relativistically moving sheath. The theoretical dis- 
persion relation allows investigation of effects associated with maximum possible 
sound speeds, Alfven wave speeds near light speed and relativistic sheath speeds. 
The prediction of increased stability of the weakly magnetized system resulting 
from c/2 sheath speeds and the stabilization of the strongly magnetized system 
resulting from c/2 sheath speeds is verified by the numerical simulation results. 
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— methods: analytical — MHD — relativity — instabilities 
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1. Introduction 

Relativistic jets arc associated with galaxies and quasars (AGN), with black hole binary 
star systems, and are thought responsible for the gamma-ray bursts (GRBs). In AGN and 
microquasar jets proper motions of intensity enhancements indicate motions that are mildly 
superluminal for the microquasar jets ~ 1.2 c (Mirabel & Rodriquez 1999), and range from 
subluminal c) to superluminal (< 6 c) along the M87 jet (Biretta et al. 1995, 1999), up 
to ~ 25 c along the 3C345 jet (Zensus et al. 1995; Steffen et al. 1995), and have inferred 
Lorentz factors 7 > 100 in the GRBs (e.g., Piran 2005). The various observed proper motions 
along AGN and microquasar jets imply speeds from ~ 0.9 c up to ~ 0.999 c, and the inferred 
speeds for the GRBs are ~ 0.99999 c. 

Jets at the larger scales may be kinetically dominated and contain relatively weak mag- 
netic fields, e.g., equipartition between magnetic and gas pressure or less, but the possibility 
of much stronger magnetic fields certainly exists closer to the acceleration and collimation 
region. Here general relativistic magnetohydrodynamic (GRMHD) simulations of jet for- 
mation (e.g., Koide et al. 1998, 1999, 2000; Nishikawa et al. 2005; De Villicrs ct al. 2003, 
2005; Hawley & Krohk 2006; McKinney & Gammie 2004; McKinney 2006; Mizuno et al. 
2006b) and earlier theoretical work (e.g., Lovelace 1976; Blandford 1976; Blandford & Zna- 
jek 1977; Blandford & Payne 1982) invoke strong magnetic fields. Additionally, Vlahakis 
and Konigl have argued that magnetically dominated "Poynting flux" jets could produce the 
accelerations observed in AGN jets such as that in NGC 6251 and 3C 345 (Vlahakis & Konigl 
2004) or provide the impetus for high Lorentz factor gamma-ray bursts outflows (Vlahakis & 
Konigl 2003). In these cases acceleration occurs up to the point at which Poynting fluxes and 
kinetic energy fluxes become comparable. In addition to strong magnetic fields, a number 
of GRMHD simulation studies of jet formation indicate that highly coUimated high speed 
jets driven by the magnetic fields threading the ergosphere may themselves reside within 
a broader wind or sheath outflow driven by the magnetic fields anchored in the accretion 
disk (e.g., McKinney 2006; Hawley & Krohk 2006; Mizuno et al. 2006b). This configura- 
tion might additionally be surrounded by a less coUimated accretion disk wind from the hot 
corona (e.g., Nishikawa et al. 2005). 

Recent observations of high speed winds in several QSO's with speeds, ~ 0.1 — 0.4c, 
also indicate that a highly relativistic jet could reside in a high speed wind or sheath, at 
least close to the central engine (Chartas et al. 2002, 2003; Pounds et al. 2003a, 2003b; 
Reeves et al. 2003). For some time other observational evidence such as limb brightening has 
been interpreted as evidence for a slower external flow surrounding a faster jet spine, e.g., 
MknSOl (Giroletti et al. 2004), M87 (Perlman et al. 2001), and a few other radio galaxy 
jets (e.g.. Swain et al. 1998; Giovannini et al. 2001). Additional circumstantial evidence 



-3- 



such as the requirement for large Lorentz factors suggested by the TeV BL Lacs when 
contrasted with much slower observed motions has been used to suggest the presence of a 
spine-sheath morphology (Ghisellini ct al. 2005). Sicmignowska et al. (2006) have proposed a 
two component (spine-sheath) model to explain the broad-band emission from the PKS 1127- 
145 jet. Additionally, a spine-sheath jet structure has been proposed based on theoretical 
arguments (e.g., Sol et al. 1989; Henri & Pelletier 1991; Laing 1996; Meier 2003) and similar 
type structure has been investigated in the context of GRB jets (e.g., Rossi et al. 2002; 
Lazzatti & Begelman 2005; Zhang et al. 2003, 2004; Morsony et al. 2006). 

In most jet generation numerical and theoretical work it has been necessary to assume an 
axisymmetric configuration. Given the helicity in the real system, e.g., helical magnetic field 
or outwards flow combined with rotation, there exists a potential problem with the stability 
of the system. Obviously, at least for most AGN jets, stability problems are surmounted and 
a highly coUimated relatively stable flow is produced. In this paper we begin a 3D numerical 
study of the stability properties of highly relativistic jet flows allowing for the effects of 
strong magnetic fields and relativistic fiow in a sheath around the highly relativistic jet. 
We note that observed relatively stable jet fiow along with observed jet structures might 
then be used to constrain the configuration in the acceleration and coUimation region where 
magnetic field strengths are high. 

In the past, 3D numerical simulations of relativistic unmagnctized jets along with ac- 
companying theoretical work (e.g., Hardee et al. 2001; Agudo et al. 2001; Hardee & Hughes 
2003) has provided relatively unambiguous interpretation and understanding of structures 
observed in the numerical simulations. Thus, we begin our study by adopting a simple 
system, no radial dependence of quantities inside the jet and also no radial dependence of 
quantities outside the jet. This "top hat" configuration can be described exactly by the 
linearized RMHD equations. In general, the system consisting of a jet with "top hat" pro- 
file and magnetic field parallel to the flow along with a uniform external medium also with 
magnetic field parallel to the fiow is more stable than a system with magnetic and fiow 
helicity. Such a system is stable to current driven (CD) modes of instability (Istomin & 
Pariev 1994, 1996; Lyubarskii 1999) but can be unstable to Kelvin-Helmholtz (KH) modes 
of instability (Hardee 2004). This approach allows us to look at the potential KH modes 
without comphcations arising from coexisting CD modes (see Baty et al. 2004) 

This paper is organized as follows. In section 2, we describe the numerical simulation 
setup and present the results of the three-dimensional RMHD simulations of spine-sheath 
relativistic jets. In section 3, we present the theoretical dispersion relation that arises from 
a normal mode analysis of the linearized RMHD equations, present relevant analytical ap- 
proximate solutions, numerically solve the dispersion relation for the simulation parameters 
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and compare to the simulation results. In section 4 we conclude. 

2. Numerical Simulations of Spine-Sheath Jets 

2.1. Numerical Method 

In order to study the long-term stability of magnetized sheath-spine relativistic jets, 
we use the 3-dimensional GRMHD code "RAISHIN" with Cartesian coordinates in special 
relativity. The method is based on a 3+1 formalism of the general relativistic conservation 
laws of particle number and energy momentum, Maxwell equations, and Ohm's law with no 
electrical resistance (ideal MHD condition) in a curved spacetime (Mizuno et al. 2006a). The 
RAISHIN code performs special relativistic calculations in Minkowski spacetime by changing 
the metric. 

In the RAISHIN code, a conservative, high-resolution shock-capturing scheme is em- 
ployed. The numerical fluxes are calculated using the Harten, Lax, & van Leer (HLL) 
approximate Riemann solver scheme. The flux-interpolated, constrained transport scheme 
is used to maintain a divergence-free magnetic field. The RAISHIN code has proven to be 
accurate to second order and has passed a number of numerical tests including highly rela- 
tivistic cases, and highly magnetized cases in both special and general relativity (Mizuno et 
al. 2006a). 

In the simulations a "preexisting" jet fiow is established across the computational do- 
main. This setup represents the case in which the jet is in equilibrium with an external 
medium far behind the leading edge Mach disk and bow shock. We allow the jet flow to be 
surrounded by a lower-density external magnetized wind medium. For all simulations, the 
ratio of densities is Pj/ Pe = 2.0, where p is the mass density in the proper frame. The jet 
flow has Uj = 0.9165c and 7 = (1 — -u^)^^/^ = 2.5. The initial magnetic fleld is assumed 
to be uniform and parallel to the jet flow. The jet is established in static total pressure 
balance with the external magnetized wind medium. Our choice of colder jet in a hotter 
wind is representative of a jet spine in a hotter sheath or cocoon as might occur as a result of 
astrophysical jet interaction with the surrounding medium. However, the specific parameters 
have been chosen for numerical and theoretical comparison convenience. 

The computational domain is 6Rj x 6Rj x 60Rj with 60 x 60 x 600 computational zones 
(10 computational zones span Rj). We impose outfiow boundary conditions on all surfaces 
except the infiow plane at 2; = 0. 

A precessional perturbation is applied at the inflow by imposing a transverse component 
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of velocity with = O.OImj. Simulations have been performed with precessional perturba- 
tions of angular frequency ujRj/uj = 0.40 (simulation A, ujI), 0.93 (simulation B, uj2) and 
2.69 (simulation C, uj3). The simulations are halted after ~ 60 light crossing times of the 
jet radius, before the perturbation has crossed the entire computational domain. 

We have performed two sets of simulations. RHD cases are weakly magnetized (see 
Table 1). The relevant sound speeds are = 0.574c and aj = 0.511c, where the sound speed 
a is given by 

Tp 



a = 



1/2 



[p + {T/T-l){p/c^)\ ' 

with r = 13/6 as the adiabatic index appropriate to a mixture of relativistically hot electrons 
and cold baryons (Synge 1957). The relevant Alfven speeds are vas = 0.0682c and vaj = 
0.064c , where the Alfven speed va is given by 



va = 



p+{T/T -l){p/c^) + B^ 



1/2 



(2) 



Therefore the Alfven speed is much smaller than sound speed. RMHD cases are strongly 
magnetized (see Table 1). The relevant sound speeds are Oe = 0.30c and aj = 0.226c. The 
relevant Alfven speeds are VAe = 0.56c and vaj = 0.45c. In this case the Alfven speeds are 
approximately twice the sound speeds. In order to investigate the effect of an external wind, 
we have performed a no wind case {ue = 0, simulation 'n') and a mildly relativistic wind 
case (up = 0.5 c, simulation 'w'). 



2.2. Numerical Results 

Figure [T] illustrates the difference in jet structure for weak magnetization with no wind 
and with an external wind. Specifically we show results from the intermediate precession 
frequency, uj2 = uRj/uj = 0.93, (cases RHDBn and RHDBw) at simulation time t = 60. 
Precession at the jet inflow plane excites the helical KH mode which is advected down the 
jet and grows. The isovolume image shows that beyond z ~ 30Rj with no wind jet flow is 
disrupted, and the magnetic field is strongly bent and distorted. Transverse 2D slices through 
the jet axis (panels lb & Ic) show a smaller density (pressure) fluctuation associated with 
the leading edge of the helix in the presence of the external wind. Transverse 2D slices 
perpendicular to the jet axis at 2; = SORj (panels Ic & If), suggest a less distorted jet in 
the presence of the external wind. Transverse velocities shown by the arrows indicate a 
circulation around the jet associated with the helical twist that is much more regular in the 
presence of the external wind. 
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(a) RHD. no wind, (0=0.93, time=60.0 




x/Rj 



Fig. 1. — Three-dimensional isovolume image {panels a, d) and two-dimensional axial {panels 
b, e) and transverse {panels c, f) slices made at z = 30Rj and simulation time t = 60 for 
the weakly magnetized cases with precession frequency u = 0.93. Panels a, b & c are 
for no wind (RHDBn) and panels d, e & f are with a wind (RHDBw). The colors show 
the logarithm of density, white lines indicate magnetic field lines {a, d), and arrows depict 
transverse velocities. 
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(a) RMHD, no wind, (0=0.93, time=60.0 
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x/Rj 

Fig. 2. — Three-dimensional isovolume image [panels a, d) and two-dimensional axial [panels 
b, e) and transverse [panels c, f) slices made at z = 30Rj and simulation time t = 60 for 
the strongly magnetized cases with precession frequency u = 0.93. Panels a, b & c are for 
no wind (RMHDBn) and panels d, e & f are with a wind (RMHDBw). The colors show 
the logarithm of density, white lines indicate magnetic field lines [a, c?)and arrows depicts 
transverse velocities. 
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Figure [2] illustrates jet structure for the strongly magnetized cases with no wind and with 
an external wind. As in Figure 1 we show results for the intermediate precession frequency, 
u;2, (cases RMHDBn and RMHDBw) at simulation time t = 60. In the no wind case the 
helical KH mode grows but more slowly than in the weakly magnetized cases shown in Figure 
1 and does not disrupt the jet inside z ~ iORj. A weakly twisted helical flow and magnetic 
structure develops. The transverse slice at z = 30Rj (panel 2c) indicates weak interaction 
between the jet and the external medium at this distance. Some circular circulation is seen. 
In the strongly magnetized case with the external wind (RMHDBw), the helical KH mode is 
damped and can just barely be seen out to z = 35Rj in the transverse axial slice (panel 2e). 
The transverse slice at z = 30Rj indicates negligible interaction between jet and external 
medium and no circular circulation (panel 2f). 

It is immediately clear from Figure 1 that presence of the wind provides a stabilizing 
influence in the weakly magnetized case. Further comparison with Figure 2 shows that the 
presence of a strong magnetic field provides a stabilizing influence and complete stabilization 
in the presence of the strongly magnetized wind. 

To investigate simulation results quantitatively, we take one-dimensional cuts through 
the computational box parallel to the z-axis at radial distances x/Rj = 0.2, 0.5, and 0.8 on 
the transverse x-axis. The results for weakly magnetized cases with/without the external 
wind are shown in Figures 3, 4, & 5 for precession frequencies, coRj/uj = 0.4, 0.93, & 2.69 
respectively. The results for the strongly magnetized cases are shown in Figures 6 & 7 for 
precession frequencies uRj/uj = 0.93, & 2.69 respectively. In the figures, and Uy velocity 
components correspond to radial Ur and azimuthal velocity components in cylindrical 
geometry. 

In the weakly magnetized cases oscillation from the growing helical KH mode is seen 
in all cases. From the plots of radial and transverse velocities, the dominant wavelengths of 
oscillation are X/Rj ~ 13, 6, & 2 for low, intermediate, and high frequency precession with 
no wind and with a measurable lengthening to X/Rj ~ 14 for low frequency precession with 
the wind. For the high-frequency case without the external wind, the azimuthal velocity 
component suggests a beat pattern with wavelength, X^eati^'^) ^ 20-Rj. From the axial 
velocities near the jet axis, we see that jet flow is disrupted at 2; ~ 25Rj, z ~ 32Rj, and 
possibly for z > 50Rj for the low, intermediate and high precession frequency no wind cases 
respectively. Jet flow is disrupted at 2; ~ 43-Rj, z ~ 35Rj and possibly for z > 50Rj for the 
low, intermediate and high frequency precession wind cases respectively. Thus, the external 
wind reduces the growth of KH instability and delays the onset of flow disruption. The 
large dips in the axial velocity near the jet surface are caused by sideways motion of the 
jet surface. Dips in the axial velocity occur more deeply inside the jet surface for lower 



-9- 



frequency precession. The presence of the external wind reduces these effects somewhat. 

Quantitative simulation results for the strongly magnetized cases with/without the ex- 
ternal wind are shown in Figures 6 and 7. In the absence of the external wind (RMHDBn 
and RMHDCn) growing oscillation from hehcal KH instabihty is seen. However, spatial 
growth is much slower than for the comparable weakly magnetized cases. Comparison with 

the weakly-magnetized cases (RHDBn and RHDCn), shows that the magnetic field reduces 
the growth rate of KH instabilities such that any disruption of collimatcd flow will lie at 
z ^ 40-Rj. In the presence of the external wind (RMHDBw and RMHDCw) the initial 
oscillation is damped. Damping is more rapid for the intermediate frequency perturbation 
than for the high frequency perturbation. Plots of radial and transverse velocities, and the 
radial magnetic field component shown in one-dimensional cuts in Figures 6 and 7 show that 
the dominant wavelengths of oscillation are A/i?j ~ 5 & 2 for the intermediate and high 
frequency perturbations, respectively with or without the external wind. There is a possible 
beat pattern in the azimuthal velocity component with A^"j(a;2) ~ 10 Rj for z < 20 Rj in 
the no wind high frequency case. The presence of the wind rcsiilts in more readily seen beat 
patterns accompanying the damped oscillations. The beat pattern is best seen in the radial 
magnetic field amplitude (Fig. 6f and Fig. 7f). At the intermediate perturbation frequency 
the beat wavelength is A^J^j(a;2) > 20 Rj. In the high-frequency case a clear beat pattern 
with \ll^^{Ld3) ~ 10 Rj is seen at 2; < 33 Rj but disappears at larger z. 

In the next section we will compare these simulation results with theoretical predictions 
for growth and damping of the helical wave mode excited by the inlet perturbations in our 
numerical simulations. 
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Fig. 3. — Radial velocity [ux)^ azimuthal velocity {uy)^ and axial velocity (7^2) along the 
one-dimensional cuts parallel to the jet axis located at xjRj = 0.2(solid hne), 0.5(dotted 
line) and 0.8(dashed line) for low-frequency precession of a weakly magnetized jet for no 
wind (RHDAn)(left panels) and with an external wind (RHDAw) (right panels). 
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Fig. 4. — Radial velocity (ux), azimuthal velocity {uy), and axial velocity {'~fUz) along the 
one-dimensional cuts parallel to the jet axis located at x/Rj = 0.2(solid line), 0.5(dotted 
line) and 0.8 (dashed line) for intermediate frequency precession of a weakly magnetized jet 
for no wind (RHDBn)(left panels) and with an external wind (RHDBw) (right panels). 
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Fig. 5. — Radial velocity (Ux), azimuthal velocity {uy), and axial velocity {'~fUz) along the 
one-dimensional cuts parallel to the jet axis located at x/Rj = 0.2(solid line), 0.5(dotted 
line) and 0.8(dashed line) for high frequency precession of a weakly magnetized jet for no 
wind (RHDCn)(left panels) and with an external wind (RHDCw) (right panels). 
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Fig. 6. — Radial velocity {u,j.), azimuthal velocity (m^,), and radial magnetic field {Bx) along 
the one-dimensional cuts parallel to the jet axis located at x/ Rj = 0.2(solid line), 0.5(dotted 
line) and 0.8(dashed line) for intermediate frequency precession of a strongly magnetized 
jet for no wind (RMHDBn)(left panels) and with a strongly magnetized external wind 
(RMHDBw) (right panels). 
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Fig. 7. — Radial velocity (m^), azimuthal velocity (m^^), and radial magnetic field (5^) 
along the one- dimensional cuts parallel to the jet axis located at x/Rj = 0.2 (solid line), 
0.5(dotted line) and 0.8(dashed line) for high frequency precession of a strongly magnetized 
jet for no wind (RMHDCn) (left panels) and with a strongly magnetized external wind 
(RMHDCw) (right panels). 
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3. Theory and Analysis 

3.1. Dispersion Relation 

Stability of a jet spine-sheath or jet-wind configuration, where the sheath or wind is 
much broader than the spine or jet, can be accomplished by modeling the jet/spine as a 
cylinder of radius R embedded in an infinite wind/sheath. 

Formally, the assumption of an infinite sheath means that the analysis could be per- 
formed in the reference frame of the sheath and numerical simulations could be performed in 
the reference frame of the sheath with results transformed to the source/observer reference 
frame. However, it is not much more difficult to derive a dispersion relation and obtain 
analytical expressions in the source/observer frame and analytical solutions to the disper- 
sion relation in the source/observer frame take on simple physically revealing forms (Hardee 
2007). Additionally, this approach lends itself to modeling the propagation and appearance 
of jet structures viewed in the source/observer frame, e.g., helical structures in the 3C120 
jet (Hardee et al. 2005). 

A dispersion relation describing the growth or damping of the normal wave modes 
associated with this system can be derived if uniform conditions are assumed within the 
jet/spine, e.g., having a uniform proper density, pj, a uniform axial magnetic field, Bj — 
Bj^z, and a uniform velocity, Uj = uj^z, and if uniform conditions are assumed within the 
external sheath/wind, e.g., having a uniform proper density, Pe, a uniform axial magnetic 
field. Be — Be z, and a uniform velocity Ug = We,z. Here the jet/spine is established in static 
total pressure balance with the external wind/sheath where the total static uniform pressure 
is = Pe + Bl/8n = P* = Pj + B]/Sn. 

The dispersion relation is obtained from the finearized ideal RMHD and Maxwell equa- 
tions, where the density, velocity, pressure and magnetic field are written as p — po + pi, 
V = u + vi (we use Vq = u for notational reasons), P = Pq + Pi, and B = Bq + Bi, where 
subscript 1 refers to a perturbation to the equihbrium quantity with subscript 0. Addition- 
ally, the Lorentz factor 7^ = (70 + 71)^ ~ 7q + 27qU ■ Vi/c^ where 71 = 7qU ■ Vi/c^. It is 
assumed that the initial equilibrium system satisfies the zero order equations. Details of the 
derivation for the fully relativistic case can be found in Hardee (2007). 

In cylindrical geometry a random perturbation pi, Vi Bi and Pi can be considered to 
consist of Fourier components of the form 

/i(r, (f), z, t) = /i(r) eiq)[i{kz ± n0 - ujt)] (3) 

where the zero order fiow is along the z-axis, and r is in the radial direction with the jet/spine 
bounded by r — R. In cylindrical geometry n is an integer azimuthal wavenumber, forn > 
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waves propagate at an angle to the flow direction, and +n and —n give wave propagation in 
the clockwise and counter-clockwise sense, respectively, when viewed in the flow direction. 
In equation (1) n = 0, 1, 2, 3, 4, etc. correspond to pinching, helical, elliptical, triangular, 
rectangular, etc. normal mode distortions of the jet, respectively. Propagation and growth 
or damping of the Fourier components can be described by a dispersion relation of the form 



XjJniPjR) XeHi'\f3,R) 



(4) 



In the dispersion relation J„ and Hn^ are Bessel and Hankel functions, and the primes denote 
derivatives of the Bessel and Hankel functions with respect to their arguments. In equation 

(4) 

— 2 2 -tjT ( 2 2 2 \ 

and 



(5a) 
(5b) 



7)2 .fTT? 



In equations (5a & 5b) and equations (6a & 6b) ro^g = {cu — kuj^ef' and g = [k 



LOU 



(6a) 



(6b) 

^e/c^)^ 



V, g/c^)"^/^ is the flow Lorentz factor, ^Aj,e = (1 - v\j^^/c^)~^/'^ is the Alfven 
Lorentz factor, W = p + [F/ (F — 1)] P/c^ is the enthalpy, a is the sound speed, is the 
Alfven wave speed, and Vms is a magnetosonic speed. The sound speed is defined by 



VP 



1/2 



where 4/3 < F < 5/3 is the adiabatic index. The Alfven wave speed defined by 



Va = 



nl/2 



1 + Vl/c^ 



where = Bq/{A7tWo) is equivalent to equation (2). A magnetosonic speed corresponding 
to the fast magnetosonic speed for propagation perpendicular to the magnetic field (e.g., 
Vlahakis & Konigl 2003) is defined by 



,1/2 
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Each normal mode n contains a single fundamental wave {cu ^ 0, k ^ 0, cu/k > 0) and 
multiple body wave{uj 0, k > 0, uj/k —>■ 0) solutions that satisfy the dispersion relation. 
In the numerical simulations the jet/spine was precessed in a manner designed to trigger 
the n = 1 fundamental helical mode. Previous theoretical and simulation work has shown 
a sufficiently close coupling between fundamental and body modes that the simulation may 
excite the first body mode as well. In what follows here we consider the n = 1 helical 
fundamental and first body mode solutions to the dispersion relation as being relevant to 
the numerical simulations performed here. 



3.2. The Helical Mode 

In the low frequency limit the helical fundamental mode has an analytic wave solution 
given by 

a; ^ [yuj + Ue] ± 17]^'^ [{uj - Uef - Vlsh^jll] ^''^ 

k (1 + vih^^c^) + vil + vy^]c^) ^ ^ 

where rj = ■jjWj /■jgWe and a "surface" Alfven speed is given by 

Vi.-bi,lF,+.Ll*'jg±g. (8) 

In equation (8) note that the Alfven Lorentz factor can be written as ^Aj,e = 1 + ^Aj,e/^'^- 
The jet is predicted to be stable to the helical fundamental mode when 

{uj - Uef - Vhh]ll < . (9) 

Thus, as might be anticipated, the growth rate is directly related to the difference between 

the magnitude of a "shear" speed, [uj — Uef, and a "surface" Alfven speed. Note that 
the "surface" Alfven speed can be greater than the speed of light and is not a physical 
wave speed. The growth rate is also reduced by the spine Lorentz factor through r] in the 
denominator of eq. (7). Finally, the real part of eq. (7) directly provides an estimate of the 
increase in helical pattern speed resulting from the external sheath fiow, and this increase 
when combined with a decrease in the temporal growth rate implies an increase in the spatial 
growth length. 

In the low frequency limit the real part of the first helical body wave solution has an 
analytic solution given approximately by 



kR ^ k'^^'R = Itt 
4 



2 2 2 2 

^]{u]-a]){u]-v%) 



1/2 

(10) 
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In this low frequency limit the body wave solution exists only when k^™R has a positive 
real part. This requires that 



'A3> 



>0. (11) 



Thus, the first helical body mode exists when the jet is supersonic and super- Alfvenic, 
i.e., m| — > and -u^ — v\- > 0, or in a limited velocity range given approximately by 
Oj > u'j > [75j/(l + 7sj)]'^i when v\j ~ a^, where = (1 — Oj/c^)'^^^ is a sonic Lorentz 
factor. 

For a supermagnetosonic jet, the helical fundamental and first body modes can have a 
distinct maximum in the growth rate at some resonant frequency. The resonance condition 
can be evaluated analytically in either the fluid limit where a ^ Va or in the magnetic limit 
where Va a. Note that in the magnetic limit, magnetic pressure balance implies that 
Bj — Be- In these cases a necessary condition for resonance is that 

1 - UjUel C2 1 + Vy,jVy,e/c^ 

where Vyjj = {aj^VAj) and v^e = {0'e-,VAe) in the fluid or magnetic limits, respectively. This 
necessary condition for resonance indicates that wc arc supersonic or super-Alfvenic when 
the shear speed exceeds a physical "surface" wave speed. When this condition is satisfled it 
can be shown that the wave speed at resonance is 

^ _ ^* ^ lj{lweVwe)Uj + -fe{lwjVwj)Ue ^^^^ 

where = (l~'i^^/c^)~^^^ is the sonic or Alfven Lorentz factor accompanying v^j = (%, vaj) 
and Vwe = ((icVAe) in the fluid or magnetic limits, respectively. The resonant wave speed 
and maximum growth rate occur at a frequency given by 

uRjv^^ ^ uj^R/vy,e = ^ • (14) 

[(1 - We/<) - {Vwe/vl - UeVy,e/c^) \ 

In equation (14) m = 0, 1 specifies the fundamental and first body modes, respectively. A 
resonant wavelength is given by A = 27ri'*/a;^ and can be calculated from 



A™ = 



m 



37r/4 + rrni 



( — ) [iyl, - Uef - [Vy,e - {Vy,eUe/ C^)vl]^\ ' R . (15) 



The resonant frequency is found to be largely a function of the sound and Alfven wave 
speeds in the sheath and the shear speed, Uj — Ug (Hardee 2007). The resonant frequency 
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increases as the sound and Alfven wave speeds increase and as the shear speed dechnes. In 
the hmit 

1 - UjUe/ 1 + VyjjV^e/ 

the resonant frequency cu^R/vwe oo. In general, When the sound or Alfven wave speed 
increases relative to the jet speed there is an increase in the growth rate at the higher 
resonant frequency accompanying an increase in the sound or Alfven wave speed relative to 
the jet speed. On the other hand, the growth rate at resonance decreases as the shear speed, 
Uj — Ue, declines. This decline in the growth rate is also indicated by equation (8) which 
applies to fundamental mode frequencies up to an order of magnitude below resonance. As 
the resonant frequency increases equation (8) applies to increasingly higher fundamental 
mode frequencies. 

Numerical solution of the dispersion relation is necessary to obtain accurate values for 
growth or damping rates as fundamental mode frequencies approach and exceed the resonant 
frequency. In general, the behavior of growth or damping associated with the first body mode 
must be obtained by numerical solution of the dispersion relation at all frequencies. In the 
high frequency limit the real part of the fundamental and first body mode solutions to the 
dispersion relation tend towards the analytic limiting form 

CU Uj ± Vwj 
— 1 — <T • 

k l±VwjUj/c^ 

which describes sound waves Vu,j — aj or Alfven waves Vy,j — vaj propagating with and 
against the jet flow inside the jet. Note that at high frequencies waves propagate in the 
spine fluid with speeds that are independent of the surrounding sheath and are decoupled 
from the spine sheath boundary. 

3.3. Numerical Solution to the Dispersion Relation 

In general, equations (7) and (10) provide initial estimates at low frequencies to the 
helical fundamental and first body mode solutions that can then be followed by root finding 
techniques to higher frequencies. The results of numerical solution to the dispersion relation 
for the parameters appropriate to the numerical simulations shown in section 2 are displayed 
in Figure 8. It should be noted that not all possible solutions are shown or have necessarily 
been found by the root finding technique. In general, in the weakly magnetized case funda- 
mental (S) mode solutions consist of a growing (shown) and damped (not shown) solution 
pair with comparable growth and damping rates (see eq. 7) and first body (Bl) mode solu- 
tions consist of a real and growing or damped solution pair. The presence of the external 
wind flow leads to reduced growth of the S mode and weak damping of the Bl mode. 
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Fig. 8. — Solutions of the dispersion relation for helical fundamental (red lines) and first 
body (green lines) modes for weakly magnetized (a^^e ^ '^Aj.e) and strongly magnetized 
(oj e ~ 'VAi,e) jet simulations without a surrounding outflow (we = 0) and with a surrounding 
0.5 c outflow [ue = 0.5). Dispersion relation solutions show the real, krRj, (dashed lines) 
and imaginary, kiRj, (dash-dot lines) parts of the dimensionless wavenumber normalized by 
the jet radius, Rj, as a function of the dimensionless angular frequency, uRj/uj, normalized 
by the jet radius and jet speed, Uj. Where the imaginary part of the wavenumber is shown 
in blue, the solution is damped. Immediately under the panel showing a dispersion relation 
solution for fundamental (S) and first body (Bl) modes is a panel that shows the wavelength, 
X/Rj, (dash-dot lines) and wave speed, v^/c, (dotted lines). The angular driving frequencies 
used in the numerical simulations are indicated by the vertical solid lines. 
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The solution structure is more complex in the strongly magnetized case. Comparable 
sound and Alfven wave speeds have led to an increased complexity compared to the weakly 
magnetized case solutions or when compared to the analytically predicted results for Alfven 
wave speed greatly exceeding the sound speed. In the absence of magnetized sheath flow S 
mode solutions again consist of a growing (shown) and damped (not shown) solution pair. 
Now however, we find multiple growing solutions associated with the Bl mode evident at the 
lower frequencies (see the lower left panel in Figure 8) . A modest damping rate accompanies 
the crossing of the real part of these multiple body mode solutions. The complex structure 
of the first body mode solution shown in Figure 8 for the strongly magnetized cases has been 
seen previously in a non-relativistic stability analysis (see Figure 20 in Hardee et al. 1995). 
At the lowest frequencies, uo < cjI, the first body mode consists of a solution pair whose real 
part can be identified with the descending and constant real part of the wavenumber. At 
frequencies above the crossing point, c<j > cjI, the rapidly rising real part of the wavenumber 
connects to the second body mode at a frequency u) > uj2. 

At the higher frequencies the Bl mode is similar to the weakly magnetized case. In the 
presence of magnetized sheath fiow, there is significant difference in growth and damping 
rates for the S mode solution pair. Weak growth is associated with the slower, Sg, moving 
shorter wavelength solution and weak damping is associated with the faster, Sf, moving 
longer wavelength solution. At the intermediate frequency, a;2, and below the growth rate is 
larger than the damping rate but at the higher frequencies somewhat above a;2 the damping 
rate is larger than growth rate for the S mode solution pair. The presence of magnetized 
sheath flow leads to damping of the Bl mode at the lower frequencies where in the absence 
of magnetized sheath flow there was modest growth and a modest high frequency damping 
rate maximum is seen where Bl intersects Ss- 

The wavelengths and growth/damping lengths normalized to the jet radius for the pre- 
cession frequencies ujl, uj2 and u3 used in the weakly magnetized simulations are given in 
Table 2. Weak damping {wd) indicates a damping length longer than the grid length in the 
numerical simulations. No entry in the growth/damping length column indicates a purely 
real solution. Wavelengths observed in the simulations are in excellent agreement with the 
theoretically predicted wavelengths. 

When no external wind is present the dispersion relation solutions show growth of the 
fundamental (S) mode at precession frequencies ujI and uj2 and growth of the body mode 
at higher frequencies. In the weakly magnetized case coupling between fundamental and 
first body (Bl) modes is indicated by the solution structure just above precession frequency, 
cuS, where the real and imaginary parts of the S and Bl mode solutions are comparable. 
Thus, we expect to see an indication of interaction in the ujS weakly magnetized simulation 
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between the S and Bl modes with a beat wavelength of Xbeat — 16.8 Rj. This agrees with 
the observed A^g„((a;3) < 20 Rj in the simulation. 

External wind flow in the weakly magnetized case leads to weak damping of the Bl 
mode and reduces the growth rate and increases the growth length, i = k~^, of the S 
mode by about a factor of two at all frequencies. However, the numerical simulation for 
the high frequency precession of a weakly magnetized jet (RHDC), see Figure 5, indicates a 
larger perturbation growth for the wind case than for the no wind case where the dispersion 
relation solutions indicate faster growth for the no wind case. The most likely reason for this 
difference is a non-linear surface and first body mode interaction in the no wind case that is 
indicated by the observed beat pattern in the no wind high frequency simulation. Here the 
different radial structure of the fundamental and body mode, see Hardee et al. (2001), with 
comparable wavelength could be responsible for destructive interference and the reduced 
transverse velocity growth seen in the no wind simulation. 

The wavelengths and growth/damping lengths normalized to the jet radius for precession 
frequencies a;l, uj2 and a;3 are given in Table 3 for the strongly magnetized simulation param- 
eters. Weak damping {wd) or weak growth {wg) indicate a damping or growth length longer 
than the grid length in the numerical simulations. No entry in the growth/damping length 
column indicates a purely real solution. Wavelengths observed in the strongly magnetized 
simulations are also in excellent agreement with the theoretically predicted wavelengths. 

When no external wind is present the dispersion relation solutions for the strongly 
magnetized case show growth of both the S and Bl modes. The growth rate of the S mode 
is reduced by about 25% at frequencies ujl and a;2 and by over a factor of three at a;3 when 
compared to the weakly magnetized case. Thus, we expect to see the observed increased 
spatial growth length in the strongly magnetized no wind simulation when compared to the 
weakly magnetized case. The Bl mode now grows at frequencies uol and uj2 in addition to 
growth at frequency a;3. Growth of the Bl mode at uil is comparable to growth of the S 
mode but is much less than the S mode at uj2 so would not be expected to appear in the 
simulation. Note that S and Bl mode wavelengths and growth rates at a;3 are comparable 
so we expect an S mode and Bl mode interaction at a beat wavelength X^eat — 10.6 Rj. This 
agrees with a weak beat pattern in and (see Figure 7) at A^"j(a;3) ~ 10 Rj seen in 
the strongly magnetized no wind a;3 simulation. 

The strongly magnetized external wind flow leads to a reduced growth rate of the 
Sg mode by over an order of magnitude at all frequencies. The damping rate of the Sf 
mode at frequencies uo < uj2 is reduced more than the Sg growth rate by about a factor 
of two. Note that for the all the other cases the low frequency growth rate of the Sg 
and damping rate (not shown) of the Sf modes are comparable. At the lower frequencies. 
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Lo < Lol the Bl mode is damped. At frequencies a;2 < a; < a;3 the Bl mode is either 
purely real or is weakly growing. At the intermediate precession frequency the strongly 
magnetized wind case numerical simulation showed damping and a beat pattern in (see 
Fi gure 6) with A^^^(u;2) ~ 20 Rj. The beat pattern can be understood as interaction 
between the weakly growing Ss and more weakly damped Sf modes with predicted beat 
wavelength \heat = 23.7 Rj. At the high precession frequency the strongly magnetized 
wind case simulation showed slower damping and a beat pattern in (see Figure 7) with 
A^^((a;3) ~ 10 Rj. The beat pattern can be understood again as interaction between the 
weakly growing Sg and now more strongly damped Sj modes with predicted beat wavelength 
\e.at — 10.0 Ry 

At the intermediate and high frequencies the dispersion relation solutions indicate weak 
growth of the Sg mode but the simulations suggest damping. However, the radial magnetic 
field component in Figures 6 & 7 show a beat pattern indicating an interaction between the 
weakly growing and weakly damped S modes. Normally, a high damping rate of the 5*/ mode 
would eliminate observable interaction. We suggest that the observed damping at uj2 and the 
lesser damping at uj?) is partially a result of this interaction. The lesser damping at a;3 occurs 
as interaction is reduced because of the considerably larger damping rate of the Sj mode at 
the higher frequency. It is also possible that some of the observed damping is a result of 
numerical dissipation given the relatively low numerical resolution of the simulations. 

We can attempt to quantify the growth or damping of the perturbations seen in the 
numerical simulations and compare the observed rates relative to theoretical predictions. 
Our estimates of the growth or damping e-folding length determined from the simulations 
are given in Table 4. The estimates are obtained by comparing a perturbation amplitude 
in determined at Z\ with a perturbation amplitude determined at z^- The e-folding 
growth or damping length I is found from A^jA^ — exp[(2;2 — -^i)/-^] ■ We always choose 2^1 > 3 
to minimize inlet effects. The range in z over which a growth or damping length was esti- 
mated is included in parentheses following the e-folding length. An indication of non-linear 
effects such as an observed beat pattern or amplitude saturation is indicated by an asterisk. 
In these cases the values provide no more than qualitative guidelines. The low frequency 
result for the no wind case is suspect as a shorter wavelength perturbation dominates growth 
at 2; > 22 and the high frequency cases all involve non-linear mode coupling, indicated by a 
beat wavelength, or appear amplitude saturated. In general, the estimates vary qualitatively 
hke the theoretical predictions. 

We can perform a more quantitative comparison of the moderate frequency, a; = a;2, 
results with theoretical growth/damping rate predictions. The simulation growth lengths 
are 1.7 — 2.5 times longer than theoretically predicted growth lengths. The weak damping 
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observed in the magnetized wind case at z > 26 when compared to the very weak growth 
predicted is consistent with this picture. Part of the difference between theory and simulation 
growth lengths can be the result of development of a shear layer in the simulations where a 
sharp boundary between spine and sheath is assumed theoretically. Additional differences 
can be a numerical effect resulting from both the narrow width of the computational domain 
and the numerical resolution. The narrow width of the computational domain means that 
some interaction with the domain boundary is unavoidable. Our relatively low numerical 
resolution of 20 computational zones across the jet diameter has the effect of increasing the 
numerical viscosity and resistivity. We do expect this to affect spatial growth lengths in 
the manner observed. We note that resolution studies in 2D RHD simulations performed 
by Perucho et al. (2004a. 2004b) suggest that only extremely high numerical resolution will 
recover correct quantitative growth or damping lengths. Excellent quantitative agreement 
between theoretically predicted and our numerically observed wavelengths indicates that we 
are not experiencing serious resolution or boundary effects. 



4. Conclusion 

We have performed numerical simulations of weakly and strongly magnetized relativistic 
jets embedded in a weakly and strongly magnetized stationary or mildly relativistic (0.5c) 
sheath using the RAISHIN code (Mizuno et al. 2006a). In the numerical simulations a jet 
with Lorentz factor 7 = 2.5 was precessed to break the initial equilibrium configuration. 
Results of the numerical simulations were compared to theoretical predictions from a normal 
mode analysis of the linearized RMHD equations describing a uniform axially magnetized 
cylindrical relativistic jet embedded in a uniform axially magnetized moving sheath. 

In the fluid limit the present simulation results confirm earlier results obtained by Hardee 
& Hughes (2003), who found that the development of sheath flow around a relativistic jet 
spine explained the partial stabilization of the jets in their numerical simulations. Here we 
confirm this earlier result and have extended the investigation to the influence of magnetic 
flelds with simulations speciflcally designed to test for stabihzation of the relativistic jet spine 
by strong magnetic fields and a weakly relativistic wind. The prediction of increased stability 
of the weakly-magnetized system with mildly relativistic sheath flow and the stabilization 
of the strongly-magnetized system with mildly relativistic sheath flow is verifled by the 
numerical simulation results. 

The simulation results show that theoretically predicted wavelengths and thus wave 
speeds are relatively accurate. On the other hand, growth rates and spatial growth lengths 
derived from the linearized equations can only be used to provide broad guidelines to the 
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rate at which perturbations grow or damp. Nevertheless, the present results can be extended 
to other parameter ranges with reduced growth occurring as 

and stabilization occurring when (eq. 9) 



In the above expressions 



represents a "surface" Alfven speed (see eqs. 7 & 8), 7Aj,e and 7j^e are Alfven and flow Lorentz 
factors, respectively; Uj^^, Bj ^ and Wj^e are the flow, axial magnetic fleld and enthalpy in the 
(j) jet or (e) external sheath. 

Formally, the present results and expressions apply only to magnetic flelds parallel to 
an axial spine-sheath flow in which conditions within the spine and within the sheath are 
independent of radius and the sheath extends to infinity. A rapid decline in perturbation 
amphtudes in the sheath as a function of radius is governed by the Hankel function's radial 
dependence. This suggests that the present analysis will provide a reasonable approximation 
to a finite sheath provided the sheath is more than about three times the spine radius in 
thickness. 

In the present regime where fiow and magnetic fields are parallel, current driven (CD) 
modes are stable (Isotomin & Pariev 1994, 1996). However, in the strong magnetic field 
regime we expect the magnetic fields in realistic jets and sheaths to have a significant toroidal 
component and an ordered helical structure. Provided radial gradients in magnetic fields 
and other jet spine/sheath properties are not too large we might expect the present results 
to remain approximately valid where uj^e and Bj^e refer to the axial or poloidal velocity 
and fleld components only. This conclusion is suggested by theoretical results, albeit non- 
relativistic and for a two dimensional slab jet, indicating that a critical parameter governing 
KH stabilization is the difference between the projection of the velocity shear and the Alfven 
speed on the normal mode wavevector (Hardee et al. 1992). In the work presented here 
magnetic and fiow field are parallel and project equally on the wavevector which for the 
hehcal mode hes at an angle 9 — tan~^{l/kR) relative to the jet axis. 

If fiow and magnetic fields are not parallel, the projection of fiow velocity and Alfven 
velocity on the wavevector is different and this will modify the stability condition somewhat. 
Of course, hehcal magnetic fields (Appl & Camenzind 1992), axially magnetized jet rota- 
tion in the subsonic limit (Bodo et al. 1996), and a radially stratified axial velocity profile 
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(Birkinshaw 1991) do modify the KH modes. Nevertheless, in the hehcally twisted magnetic 
and flow field regime likely to be relevant to most astrophysical jets where CD modes are 
unstable (Lyubarskii 1999), there can be competition between CD and KH modes. At least 
in the force-free magnetic field regime, KH modes can dominate CD modes when both are 
unstable (Appl 1996). 

While the normal Fourier modes, such as the helical mode that we have considered in 
this work, are the same in KH and CD regimes, the conditions for instability, the radial 
structure, the growth rate and mode motions are different. Non-relativistic simulation work 
(e.g., Lery et al. 2000; Baty & Keppens 2003; Nakamura & Meier 2004) suggests that CD 
structure is internal and moves with nearly the jet speed. On the other hand, KH structure 
is surface driven and can move at speeds much less than the jet speed. These differences 
may serve to identify the source of helical or other moving structure on relativistic jet flows 
and allow determination of jet properties near to the central engine required to produce such 
structure. 
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Table 1. Models and Parameters 



Case 


ujRj/uj 




Oe/C 


aj/c 


VAe/c 


VAj/c 


RHDAn 


0.40 


0.0 


0.574 


0.511 


0.0682 


0.064 


RHDBn 


0.93 


0.0 


0.574 


0.511 


0.0682 


0.064 


RHDCn 


2.69 


0.0 


0.574 


0.511 


0.0682 


0.064 


RHDAw 


0.40 


0.5 


0.574 


0.511 


0.0682 


0.064 


RHDBw 


0.93 


0.5 


0.574 


0.511 


0.0682 


0.064 


RHDCw 


2.69 


0.5 


0.574 


0.511 


0.0682 


0.064 


RMHDBn 


0.93 


0.0 


0.30 


0.226 


0.56 


0.45 


RMHDCn 


2.69 


0.0 


0.30 


0.226 


0.56 


0.45 


RMHDBw 


0.93 


0.5 


0.30 


0.226 


0.56 


0.45 


RMHDCw 


2.69 


0.5 


0.30 


0.226 


0.56 


0.45 



Table 2. Wave & Growth/Damping lengths: Weakly Magnetized 



LoRj/Uj 


A5(0) 


^5(0) 


Ab(0) 


MO) 






><B{Ue) 


^B{Ue) 


0.40 


13.0 


6.85 


4.3 




14.5 


U.7g 


4.4 


wd 


0.93 


5.4 


3.35 


3.2 




6.2 


6.9(7 


3.1 


wd 


2.69 


1.83t 


1.9g 


1.65^ 


3.3g 


2.14 


2.7g 


1.58 


wd 



^(^beat) ^ - ^B^ - 



Table 3. Wave & Growth/Damping lengths: Strongly Magnetized 



uRj/uj A5(0) £5(0) Ab(0) ^b(O) As(Ue) ^sK) \B{Ue) isiUe) 



0.40 



0.93 



12.8 9.75 



5.4 5.75 



4.65 
5.04 



2.38 
3.15 



2.69 1.94t 6.25 1.64t 



10.95 
7Ag 



11.7 
15.1 



wg 
wd 





4.2 




38d 




4.6 




wd 



5.1 
6.5 



625 
wd 



1.77 




wg 


2.15 




19.6d 



3.1 



1.69 



Table 4. Simulation Growth/Damping e-folding lengths 



ujRj/uj 


RHDn 


RHDw 


RMHDn 


RMHDw 


0.40 


26 (4-22) 


30 (3-40) 






0.93 


9 (5-22) 


12 (4-38) 


14 (5-37) 


> 25* wd (26-43) 


2.69 


14* (11-27) 


23* (7-27) 


> 25* wg (11-25) 


> 19* wd (18-35) 



*mode interaction or saturation 



